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Abstract 

Response construction methods using Moving Least Squares (MLS), Kriging and Radial Basis Functions (RBF) are 
compared with the Global Least Squares (GLS) method in three numerical examples for derivative generation 
capability. Also, a new Interpolating Moving Least Squares (IMLS) method adopted from the meshless method is 
presented. It is found that the response surface construction methods using the Kriging and RBF interpolation 
yields more accurate results compared with MLS and GLS methods. Several computational aspects of the response 
surface construction methods also discussed. 


I. Introduction: 

Structural reliability engineering analysis involves determination of the probability of structural failure taking into 
account the uncertainty in the geometric parameters, material properties and loading conditions [1]. The uncertain 
quantities are treated as random variables, and often a large number of simulations (structural analyses) with 
different sets of random variables are necessary to estimate the reliability of a structure. Hence, the computational 
effort required to perform a structural reliability analysis can be very high. In order to minimize the computational 
time, response surface functions are often used as simple and inexpensive replacements for computationally 
expensive structural analyses in reliability methods. Most of the response surface construction methods use Global 
Least Squares (GLS) methods. The GLS methods use a single quadratic or cubic polynomial representing the entire 
parametric space of the random variables. In general, a single polynomial to represent the entire parametric space 
introduces large errors in the response estimation or limits the size of the parametric space. In order to overcome 
such difficulties, local (piecewise) polynomial functions with higher order derivative continuity were used in 
response surface approximation [2-4]. In arriving at an interpolated value at some point in the parametric space, the 
local methods more heavily weight data samples that are “nearby” rather than giving all data samples equal weight. 
It is found that higher order derivative continuity methods require fewer sampling points for the same accuracy [2] 
compared to GLS methods. 

A new local Moving Least Squares (MLS) response surface construction method was developed in reference 2. The 
MLS method was compared with other local methods such as Kriging [3] and found to be more accurate and 
computationally effective for the examples considered. Another class of local response surface construction using 
Radial Basis Functions (RBF) was proposed in reference 4. Along with the classical RBF, polynomial augmented 
RBF and compactly supported radial basis functions [5-7] were used for response surface construction in reference 
4. The methods based on MLS, Kriging, and RBF were compared in two numerical examples in references 2 and 4 
and found to predict the response almost with the same accuracy. Even though the MLS, Kriging, and RBF methods 
produce higher accuracy in response prediction compared to the GLS methods, the computational efficiency or the 
ability of these local methods to reproduce derivatives is not investigated in the literature. The ability to estimate the 
derivatives accurately is very critical for the response construction methods in optimization and reliability studies. 
Hence, it is necessary to examine the local methods for derivative generation capabilities. The purpose this paper is 
to study and compare the computational efficiency and derivative generation capability of the local MLS, Kriging 
and RBF methods by applying these methods to several numerical examples. Also, a new Interpolating Moving 
Least Squares (IMLS) method adopted from the meshless method is presented. 

First, a brief introduction to MLS, RBF, Kriging and GLS methods of response surface construction is presented. 
Equations are presented to estimate the first partial derivative for all the methods. Second, the selected response 
surface construction methods are tested in three numerical examples to demonstrate the effectiveness of the local 
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methods for response prediction and derivative estimation. Also various computational aspects of the response 
surface construction methods are discussed. 

II. Interpolation methods used for Response Surface Construction 

In this section, the interpolation methods used in the study are presented. It is required to construct a response 
surface from the given function values F i (i = 1,2, ...IV ) at N sampling point locations x,- (i = 1,2,...N ). Here x 

denotes the vector of space coordinates x=jx 1 ,x 2 ,...x rf j with d as the dimensions of the interpolating surface. 
The order of the smoothness of the generated surface is measured by the ability of the surface to produce the partial 
derivatives. The surface is C 1 continuous, if the function values and first derivatives are continuous. The surface is 

C~ continuous if the function values and its second derivatives are continuous. The surface is C continuous, if 
the function values and third derivatives are continuous everywhere in the parametric space. 

A. Moving Least Squares (MLS) method: 

The Moving Least Squares (MLS) method is widely used in meshless methods [8], The MLS approximation for 
the estimated value s(x) at an arbitrary point can be written as 


s{x}= p T {x}a m {x) 


(1) 


where p T {x} = \p\{x},p 2 {x},...p m {x}\ is a polynomial basis function of order m used in the MLS interpolation. 
The coefficients a j(x), j = l,2,...,m , are functions of the spatial coordinates. For example, for two design 
variables, 


P 


T 



Linear basis function; m = 3 


p T {x}= l,x 1 ,x 2 (x 1 ) 2 ,x 1 x 2 ,(x 2 ) 2 


Quadratic basis function; m = 6 


(2) 


The unknown coefficients a m can be determined using the weighted least squares error norm J(x) at the N 
sampling points 


J{x)='Lw i {x\p T (x i )a m {x)-Y 

i=l (3) 

= [Pa m (x)-YY -W-[Pa m (x)-Y] 

where w t (x) is the weight function associated with node i , whose value is nonzero only in the support or influence 
domain of node x, (usually a sphere of radius I , ). The matrices P and W are defined as 


P = 


P T (*i) 
P T (x 2 ) 


P T (x n ) 


Nxm 


( 4 ) 
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0 ••• w N (x) 


-I NxN diagonal matrix 


and the vector Y = \}’\,y 2 > where are the fictitious values at the N sampling points. 

The fictitious values are unknown at this point in the interpolation scheme and will be selected later. 

The unknown constants a m (x ) can be determined by minimizing the norm J(x) in Equation (3) with respect to 
a m (x) which leads to the following linear relation between a m {x) and Y , 

A(x)a m (x) = B(x)Y (6) 

where the matrices A(x) is a square matrix of m rows and m columns, and B(x) is a matrix of m rows and N 
columns. These matrices are defined by 


A(x) = P T WP = B(x)P = ^w i {x)p{x i )p T (x i ) (7) 

i=l 

B(x) = P T W = [wj (x)p(x l ), w 2 (x)p(x 2 w N (x)p(x N )] (8) 

The MLS approximation is defined only when the matrix A(x) in Equation (7) is non-singular. The matrix A(x) 
will be non-singular if at least m weight functions are non-zero within the influence domain of the node x,- . 

Substituting Equation (6) in the Equation (1), the estimated function values at an arbitrary point x can be obtained 



Figure 1. Classical Moving Least Squares (CMLS) method. 

in terms of the not yet selected fictitious values Y as 

s(x) = p T (x)A 1 (x)7?(x)T (9) 

The influence domain radius /, is a free parameter in the MLS method and pre-selected by the user. The value 
of the radius /,• affects the computational efficiency and accuracy of the MLS method. Depending upon the choice 
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of the fictitious values for , the MLS interpolation can be termed as Classical Moving Least squares (CMLS) or 
Interpolating Moving Least Squares (IMLS) as described below: 

Classical Moving Least Squares (CMLS): 

Given the function values F t (for i = 1,2,... N) at the N sampling points, the fictitious values are selected in 
CMLS such that 


y, = F t for i = 1,2, ...IV (10) 

The fitted surface in CMLS is a true least squares surface and does not reproduce the function values F i at the 
sampling points i.e., 


s(Xj ) ^ Fj for i = \,2,...N (11) 

The CMLS interpolation is shown schematically in Figure 1. In CMLS for every new interpolation point, the 
matrices in Equation (9) are evaluated. This equation involves assembly and inversion of the (in x in) matrix A and 
also assembly of the matrix B . 


Interpolating Moving Least Squares (IMLS): 

Given the function values 7y( for i = 1,2,... /V) at the N sampling points, the fictitious values are evaluated such 
that 


s(Xj ) = Fj for i = 1,2 , ...N 

Using Equation (12), the fictitious values y, are evaluated using Equation (9) as follows 

M = HY 


( 12 ) 


(13) 


where 



s ( x l ) 


1 

A 1 ( x)B( x ) 

x=x x 


y i 

M = < 

s(x 2 ) 

; H = 

P T ( *2 )\ 

A 1 ( x)B( x ) 

x=x 2 

Y = < 

T2 


s(x N ) 

(Nxl) 

p r ( x n) 

A X (x )B( x) 

X=X N _ 

(NxN) 

Yn. 


(14) 


(tVxl 


or 


Y =H~ X M 


(15) 
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Once the fictitious values are evaluated using Equation (15), the function values at an arbitrary point X can be 
estimated using Equation (9). The fitted surface in IMLS is a true interpolating surface, implying that the surface 
passes through each and every sampling point function value. The IMLS interpolation is shown schematically in 



Figure 2. Interpolating Least Squares (IMLS) method. 


Figure 2. The fictitious values obtained in the IMLS have no physical significance, except to make the interpolation 
pass through all of the sampling points. 

The IMLS method involves inverting a large matrix H of (N xN) once in the interpolation to evaluate the 
fictitious values >7 . Once the fictitious values are evaluated, the rest of the computational procedure is the same as in 
the CMLS method for interpolating the function value at an arbitrary point x . 

Selection of the weight function: 

As already mentioned, the weight function w,(x) is non-zero only in the influence domain of node i , and it is 
equal to zero outside the influence domain. The weight function is selected such that its value ranges from unity at 
the center of the influence domain to zero at the boundary and outside the influence domain. In the present study, the 
influence domain is assumed to be a sphere with radius l i . The radius /,- must be large enough to contain at least 
m nodes in each direction of the parametric space. 

12 3 

In this paper, three spline functions with C ,C , and C continuity are used as the weight function 


For C 1 : 


w’i(x) = 


l-3t[+2tf 


0 < ti < 1 
tj >1 


For C 2 


Wi(x) = 


-10t 2 


+ \5tf-6 tf 
0 


0<t t < 1 

tj >1 


(16) 


(17) 
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0<h< 1 

U > 1 


(18) 


For C 3 


Wi(x) = 


1-3517+84 tf 


-70 if +201, 6 7 


0 


pc - x t 

where 1 ; - = is the normalized distance, from the center of the influence domain ( x,-) and a general point x . 

h 

The smoothness of MLS approximation is controlled by both the weight and the basis functions. The precision 
(continuity) of MLS interpolation will be equal to the minimum precision of the weight and basis functions. 

Estimating Derivatives in MLS: 

The partial derivative estimation for the MLS method is complicated, since the coefficient a{x) in Equation (1) 
is a function of x instead of constant as in all other interpolating methods. Hence, the derivative estimation in the 
MLS method is computationally more intensive than all the methods presented here. 

The first partial derivatives of the function ,s(x) for both CMLS and IMLS can be obtained by differentiating 
Equation (9) as 


dS(x) 

dx 


dp T (x) 


dx 


\A~\x)B(x ) 


+ P T (x) I 


A -i (x) ?m + ^i B(x ; 


dx 


dx 


where 


dA~\x) 

dx 


dx 


A-\x) 


(19) 


(20) 


B. Radial Basis Function (RBF) method 

In classical RBF based methods, the interpolation of a surface s(x) is performed as a linear combination of 
radial functions as 


s(x)='£A i ip (Jx-x^lJ (21) 

1=1 

where the radial basis functions cp are functions of the radial distance ||x - x,| 2 from node i , are interpolation 
constants to be determined, and N is the number of sample or data points with known function values F i such that 

s(Xi) = Ff ( 22 ) 

The Euclidean norm |x — x,-| 2 represents the radial distance r of the point x from the center Xj . For two- 
dimensional systems in Cartesian X-Y coordinates, the radial distance can be obtained as 

|x-x,.|| 2 =j(X-X i ) 2 +(Y-Y i ) 2 (23) 
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The unknown interpolation coefficients Xj in Equation (21) can be determined by minimizing the norm 


J = 



z Ai<p (h-*4 


) 


for = 1,2,3 


( 24 ) 


The minimization equation in matrix form can be written as 

MW=W 


where 


{a} T - {X x ,A 2 ,...,X n } 

{Rf ={F h F 2 ,...,F N ) 

and the coefficient ( i th row and j th column) of the matrix [a] can be obtained from 

a ij =<p\ x i -*/|| 2 ) 


(25) 


(26) 


(27) 


Since there are N equations with N unknown constants X,- , i = 1,2 ,- --N , the resulting surface is an interpolating 
surface. The most commonly used classical radial functions are given in Table 1: 


Table L Commonly Used Classical Radial Basis Function s 


Name 

RBF Function 

r= \\ x ~ x i 2 

Linear 

<f S(r) = cr 

Cubic 

</>(r) = (r + c) 3 

Thin plate spline 

<, i>{r ) = r log(cr“) 

Gaussian 

r V 

0 

1 

II 

Multiquadratic 

(/>{>■) = ( r 2 +c 2 ) 2 


Recently, compactly supported radial basis functions are generated [5-8]. In compactly supported radial functions, 
the interpolation matrix [/)] in Equation (25) is sparse and positive definite. Also, for interpolation, only a few 
terms need to be considered making the algorithm more efficient for the computation and evaluation of response 
surfaces. In this study, the following two compactly supported radial functions were used. These functions were 
adapted from reference 5. 
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Compact-I: 


<p(r) = (1 - 0 5 (8 + 40 1 + 48 r + 25t 3 + 5 1 4 ), 

0 < r < r 0 


(28) 


= 0 


r > r 0 


Compact-II: 


<p(r) = (1 - 0 6 ( 6 + 36 1 + 82 1 2 + 72 1 3 + 30 1 4 + 5 1 5 ), 

0 < r < r 0 


(29) 


= 0 


r > r 0 


r 

where t = — , and r 0 is the radius of the domain of compact support. The compact support radius r 0 is a free 
r 0 

parameter in the interpolation, and is selected by the user. Many other compact support functions are given in 
references 5 and 6. The classical radial function given in Table 1 and the compactly supported radial basis functions 
in Equations (28) and (29) are used in Equation (21) for response surface generation. 

Augmented RBF Response Surface Construction 

The classical radial functions given in Table 1 and the compactly supported RBF in Equations (28) and (29) do 
not reproduce constant, linear polynomial terms [4], The inability to reproduce the polynomial terms is not desirable 
if the response is in terms of the simple polynomials. In order for the RBF to reproduce simple polynomials, the 
functions can be augmented by polynomial terms as 


N , m 

s(X)= XA ( #)+ I P/xjbj (30) 

i= 1 7=1 

where Pj(x) are the monomial terms in the polynomials and bj are the additional m constants introduced in the 
interpolation due to the polynomial terms. For two-dimensional problems, the monomial terms are 

n 1 2 , K2 12, 2,2 s 

(1,* ,X ,(x ) ,x X ,(x ) ,...). 

The additional m unknown constants introduced in Equation (30) depend on the order of the polynomial selected 
and the order of the problem dimension d or number of design variables. The number of m unknown constants due 
to the polynomial terms can be determined from 


m = 1 


for constant 


m = d + 1 


( d + 1)(^/ +2) 


for linear polynomial 


for quadratic polynomial 


( 31 ) 


m = ( d + l >(d + 2)(d + 3) for cub . c polynomial 
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Equation (30) consists of N equations with N + m unknowns. The additional m equations necessary to solve for 
the m additional unknowns can be obtained from the following m constraints [4] 


N 


'LPj(x)^ 


for j = 1,2,3,- 


(32) 


Hence, the polynomial augmented radial functions do not require any additional sampling points. The constraints in 
Equation (32) are handled internally in the interpolation. Equation (25) (to obtain the unknown constants) can be 
modified (using Equations (30) and (32)) for the augmented RBF as 


MyVxtV l^lvxw JWatxI I _ [{MtVxl j 

[Bl xN [oL xm J l Mm xl J [Mm xl J 


where 


(33) 


{bf ={b h b 2 ,b 3 ,...,b m } 


(34) 


and the coefficient bjj ( i th row and j th 


column) of the matrix [b\ can be obtained from 


kj = Pj(*i) 


j = 1,2, and i = 1,2,---,N 


(35) 


Since, there are N + m equations with N + m unknown constants, the resulting surface constructed using 
augmented RBF is also an interpolating surface. It is important to point out that, when augmented RBF is used to 
represent simple polynomials, the constant Aj in Equation (30) will be identically zero. For non-polynomial 

functions, the polynomial terms are expected to augment the performance of the RBF. In the present study, the RBF 
are augmented with complete cubic polynomials. 

Estimating Derivatives in RBF: 

The partial derivatives of RBF are much simpler than the MLS method, since the coefficients are constant. The first 
order partial derivatives can be obtained by differentiating Equation (30) as 


s(x) R 8</>(r) dr ™ bPj(x) 

— A A; h ?, O : 

dx i = i dr dx J= \ dx J 


(36) 


C. Kriging method 

Kriging is an interpolation method that originated in geostatistics. Kriging uses the properties of the spatial 
correlation among the data samples. In arriving at an interpolated value at some point in the parameter space, 
Kriging more heavily weights data samples that are “nearby” rather than giving all data samples equal weight. In 
Kriging. interpolation is achieved by setting the mean residual error to zero and also by minimizing the variance of 
the errors. The final equations for Kriging are given below from reference 3 for N sampling points and d design 
variables. The estimated value of s in Kriging is obtained from 

s(x) = P + r T R- l ^-Pf) (37) 
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where Y is the column vector of known function values at the N sampling points, P is a constant to be determined, 
R is a correlation matrix obtained for an i ,h row and j th column from the correlation function as 


R(x i,x j) = exp 


6 Z xf - x) 
k = 1 


(38) 


r is the column vector of length N obtained from 

r = ),R(x,x ^ ),. . .,R(x,Xpj )f (39) 

and / is vector of length N , with all elements in the vector set to unity as 

/Mu, i ,i F ( 40 ) 

The unknown /? in Equation (7) can be obtained from 

P = (f T R-\n~ l f T R- l Y (41) 

The unknown quantity 0 in Equation (38) is obtained from the Maximum Likelihood Estimate (MLE) in a one- 
dimensional maximization problem defined by 


MLE= max 


-][/V ln((j 2 ) + ln|7?| 


0 < 0 < 00 


(42) 


where 


(Y-fff R-'(Y-ff) 
N 


( 43 ) 
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Estimation of 6 in the one-dimensional optimization problem is the critical and most computationally 
expensive step in the Kriging method as it involves using an optimization code. The Kriging method used in this 

study produces a C 2 continuous interpolating function over the entire parametric space. Once 0 is known, the P 
in Equation (37) can be uniquely determined to complete the interpolation. 

Estimating Derivatives in Kriging: 

The first order partial derivatives for the Kriging surface can be obtained by differentiating Equation (37) as 

^ = (44) 

OX ox 


with 


and using 


dr T (x) _ jdR(x,Xi ) dR(x,Xi) 
dx I dx dx 


3R(x,xn ) 
dx } 


dR(x,x j ) 
dx 


d , 

= -eH{x k 

k=l 




d , 

e z ( x k 

k = 1 


ks2 

X.) 


fory = l,2,"-,7V 


(45) 


(46) 


D. Global Least Squares (GLS) method 

The GLS methods are generally known as polynomial regression methods and are the most commonly used 
method in the literature [3], The GLS approximation for the estimated value s(x) at an arbitrary point using a 
quadratic polynomial can be written as 


d d . 

s(x) = a 0 + Z « jX J + Z bjjX 1 x 1 (47) 

7=1 


where d is the dimension of the problem or number of design variables. The unknown coefficients aQ,aj, and by 

are determined by a regression procedure. Most commonly, the method of least squares is used to determine the 
coefficients by minimizing the error of the approximation at the sampling points. Since a single polynomial is used 
to represent the entire parametric space, the method is called the Global Least Squares (GLS) method. In the present 
study, the GLS method is limited to the quadratic polynomial given by the Equation (47). 

Estimating Derivatives in GLS: 

The first order partial derivatives can be obtained by differentiating Equation (47) as 


0.?(x) 
dx i 


~ a j + 


d 

z 


byX 1 


(48) 
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ITT. Application Problems: 

The five response surface generation methods described in the previous sections were tested in three numerical 
examples. The response surfaces were generated to fit the three one-dimensional curves given in Table 2 using the 
function values given at the sampling points. The number of sampling points was selected to be the minimum 
number needed to reproduce the function F(x ) with reasonable accuracy by any one of the methods. Then, the 
selected number of sampling points was used in all other methods to compare to one another. The basis function in 
the IMLS and CMLS methods are restricted to a quadratic polynomial to be consistent with the quadratic 
polynomial used in the GLS method in these problems. 


Table 2. One-Dimensional Application Problems Studied 



Equation 

First Derivative 

Parametric Space 

Number of 
Sampling 
Points Used 

Quadratic Polynomial 

F(x) = 1 + x + x" 

dF(x) _ j + 2 ^ 
dx 

0 < x < 7 

5 

Trigonometric Function 

F(x) = x cos(x) 

dF(x) 

= cox(x) - x sin(x) 

dx 

0 < x < 7 

9 

Rational Function 

x 2 

F(x) = —^ 
8 + x 5 

dF(x) 16x-3x 6 
dx (8 + x 5 ) 2 

0 < x < 7 

21 


The response surfaces generated by the different methods were evaluated by assessing the method’s ability to 
represent the functions accurately and also its ability to reproduce the first partial derivative of the design variable. 
The function values and first order derivatives were interpolated at 201 equally-spaced points in the parametric 
space. The interpolated values obtained using response surfaces were compared with the theoretical/exact values that 
were calculated using the equations given in Table 2. An average percent of error defined by 

X |(exact)y - (predicted) 

Average percent Error = — — x 100.0 (49) 

X {exact) 
i = 1 

is used to compare the accuracy of the methods. In equation (49) the “exact” values are calculated using the 
equations in Table 2, and the “predicted” values are obtained by interpolating the response surfaces. 
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A. Quadratic Polynomial Function 

In this study, a quadratic polynomial is used as the basis function in the GLS, CMLS and IMLS methods. Hence 
these methods are expected to reproduce up to a quadratic polynomial exactly. Also, the cubic polynomial 
augmented RBF functions can reproduce the quadratic polynomial exactly. It is known that the Kriging method can 

produce aC continuous surface and hence can also reproduce quadratic functions exactly. It was verified in the 
numerical examples tested that all five response surface methods reproduced the quadratic polynomial function and 
its first derivative exactly everywhere in the parametric space to the computer machine accuracy. 


B. Trigonometric Function 

The shape of the trigonometric function and the position of the nine Progressive Lattice Sampling (PLS) points 
[2] used are shown in Figure 3. In the PLS methods as the number of sampling points increased from one level to 
another, all the sampling points from the previous level are retained in the current level. 



Figure 3. Trigonometric function and position of nine sampling points. 


The estimated function values at the nine sampling points are used to generate the response surfaces for all five 
methods. Then, the responses are estimated at 20 1 equally-spaced points by interpolating the response surface and 
the estimated values are compared with the exact values. 


IMLS and CMLS Interpolation Methods for Trigonometric Function 

The responses obtained by the IMLS and CMLS methods are given in the Figures 4 and 5. For clarity in the 
presentation of the results, only 40 of the 201 interpolated values are shown in figures, herein. As seen in Figures 4 
and 5, the IMLS interpolation predicts the response and the derivative more accurately than the CMLS method. In 
obtaining the responses in Figures 4 and 5, the free parameter /,■ was selected to yield accurate values. 

It can be noted that in the IMLS and CMLS methods, it is possible to increase the order of the basis function 
very easily without large computational penalty. For the numerical examples presented in this study, for MLS 
methods, the order of the basis function can be changed to cubic polynomial to reduce the error in the interpolation 
of the function and in the estimation of derivatives. In the RBF and Kriging methods, the order of the basis function 
is fixed and can not be changed. 
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(b) 


Figure 4. IMLS interpolation of the trigonometric function: (a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative response 
values. 
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(a) 



(b) 

Figure 5. CMLS interpolation of the Trigonometric Function: (a) Comparison of function 
values and interpolated response values (b) Comparison of derivative values and interpolated 
derivative response values. 
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RBF Interpolation for Trigonometric Function 

All the classical Radial Basis Functions in Table 1 and the compact RBF in Equations (28) and (29) were used in 
the response surface construction. The free parameter c in the RBF shown in Table 1 was selected to be equal to 
unity in this study. There is no other free parameter to be selected in the RBF methods. Flowever, there is a need to 
select best type of RBF from the large class of available radial functions. The average errors obtained in 
interpolating the function values and derivative estimation using Equation (49) are shown in Table 3 for the 
functions studied herein. 


Table 3. Average error obtained for various RBF with trigonometric functions 


RBF 

Average Function Error (%) 

Average Percent Error in 
Derivative Estimation 

Linear 

7.8 

31.66 

Cubic 

1.06 

4.61 

Thin Plate 

2.25 

10.0 

Gaussian 

0.06 

0.274 

Multiquadratic 

0.1041 

0.462 

Compact I 

0.83 

3.65 

Comapct II 

0.363 

1.6 


As observed in Table 3 all the Radial Basis Functions reproduced the function values accurately within five 
percent, except for the linear function. The function values and derivatives obtained using Gaussian and 
multiquadratic radial basis functions are most accurate. The least average percentage error is obtained for the 
function values for augmented Gaussian RBF. The comparison of responses for the function and derivatives are 
shown for the Gaussian function in Figure 6. 


Kriging Interpolation for Trigonometric Function 

The response obtained using Kriging interpolation is compared in Figure 7. The free parameter 9 obtained using 
optimization has a value of 2.0 . As illustrated in Figure 7, the Kriging method interpolates the trigonometric 
function very accurately. 
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(a) 


7 


RBF Interpolation 



-5 


(b) 

Figure 6. RBF interpolation of the trigonometric function:(a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative response. 
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7 

- X Kriging Interpolation 



x 


-5 


Figure 7. Kriging interpolation of the trigonometric function: (a) Comparison of function values 
and interpolated response values (b) Comparison of derivative values and interpolated derivative 
response values. 
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GLS Interpolation Method for Trigonometric Function 

The response obtained using the GLS interpolation is shown in Figure 8. Since the GLS interpolation uses a 
single polynomial for the entire parametric space, the error obtained is large compared to all the other methods in the 
study. 
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40 



(b) 


Figure 8. GLS interpolation of the trigonometric function: (a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative response 
values. 
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Comparison of response surface methods for the Trigonometric function. 


All five methods used in the response surface construction for the trigonometric function are compared in Table 
4. The average error calculated using Equation (49) is compared in the table. 

Table 4 : Comparison of errors obtained in response surface methods for the Trigonometric function 


Response Surface 
method 

Average Function Error (%) 

Average Percent Error 
in Derivative Estimation 

IMLS 

1.2092 

6.65 

CMLS 

15.12 

26.12 

RBF 

0.060 

0.274 

Kriging 

0.0 

0.238 

GLS 

49.12 

93.86 


The Kriging and RBF methods are the most accurate surface construction methods for this problem. In the 
moving least squares, the IMLS method performs better than the CMLS. The GLS method continues to yield poor 
accuracy compared to all other methods. 

C. Rational Function 

Initial studies indicated that the minimum number of sampling points needed to reproduce the rational function is 
considerably higher than the trigonometric function in the previous section. Hence, twenty one equally-spaced 
sampling points were used. The rational function along with 21 equally-spaced sampling points is shown in Figure 
9. 
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IMLS and CMLS Interpolation Methods for Rational Function 

The responses obtained using the IMLS and the CMLS interpolation functions are shown in Figures 10 and 11, 
respectively. 



(a) 



(b) 

Figure 10. IMLS interpolation of the rational function: (a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative 
response values. 
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(a) 



(b) 

Figure 11. CMLS interpolation of the rational function: (a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative response 
values. 
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RBF Interpolation Methods for Rational Function 


All the classical RBFs in Table 1 and compact RBFs in Equations (28) and (29) were used in the response 
surface construction. The average errors obtained in interpolating the function values and derivative estimation using 
Equation (49) are shown in Table 5 for all the RBF used. As seen in Table 5, all the RBFs reproduced the function 
values accurately except the linear function. The Thin plate and Gaussian radial basis functions were very accurate. 


Table 5. Average Error Obtained using RBF with Rational Functions 


RBF 

Average Function Error (%) 

Average Percent Error in 
Derivative Estimation 

Linear 

2.37 

20.76 

Cubic 

0.67 

12.57 

Thin Plate 

0.534 

5.3 

Gaussian 

0.89 

6.1 

Multiquadratic 

1.78 

7.9 

Compact I 

1.26 

5.8 

Comapct II 

2.66 

16.0 


The RBF interpolation using with the thin plate rational function is shown in Figure 12. 
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(b) 

Figure 12. RBF interpolation of the rational function: (a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative response 
values. 
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Kriging Interpolation for the Rational Function 

The Kriging interpolation for the rational function is shown in Figure 13. The free parameter 6 obtained using 
optimization has a value of 16.0 . Here again, the Kriging method reproduced both the function and derivative for 
the rational function extremely well. 



(a) 



(b) 

Figure 13. Kriging interpolation of the rational function: (a) Comparison of function values and 
interpolated response values (b) Comparison of derivative values and interpolated derivative 
response values 
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GLS Interpolation Method for Rational Function 

The response obtained using GLS interpolation is shown in Figure 14. Since GLS uses a single polynomial for 
the entire parametric space, the error obtained is large compared to all the other methods in the study. 



(a) 



(b) 

Figure 14. GLS interpolation of the rational function: (a) Comparison of function values and interpolated 
response values (b) Comparison of derivative values and interpolated derivative response values. 
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Comparison of response surface methods for rational function. 


All the five methods used in the response surface construction for the rational function are compared in Table 6. 
The average error calculated using Equation (49) is compared in the table. The Kriging and RBF methods are the 
most accurate surface construction methods for this problem. In the moving least squares, the IMLS method 
performs better than the CMLS. The GLS method continues to yield poor accuracy compared to all other methods. 


Table 6. Comparison of Errors in Response Surface Methods for the Rational Function 


Response Surface 
method 

Average Function Error (%) 

Average Percent Error in 
Derivative Estimation 

IMLS 

0.68 

5.53 

CMLS 

1.14 

6,13 

RBF 

0.53 

5.30 

Kriging 

0.008 

0.798 

GLS 

66.95 

100.0 


IV. Discussion 

Based on the results presented in the previous section the following observation can be made 


i. The Kriging and the Radial Basis Function (RBF) methods are the most accurate methods for the 
numerical examples studied. 

ii. The Global Least Squares (GLS) method consistently behaved poorly both in function interpolation and 
in derivative estimation. 

iii. In the Moving Least Squares (MLS) methods, the Interpolating Moving Least Squares (IMLS) method 
is more accurate than the Classical Moving Least Squares (CMLS) method. 

iv. Changing the order of the basis functional (example: from quadratic to cubic) is possible only in the 
case of MLS methods. In all the other methods, the basis function can not be changed. 

v. Estimation of free parameter in the Kriging using the optimization procedure is the most expensive step. 

vi. Even though there is no free parameter needed to be selected in RBF, selecting the correct type of radial 
basis function for a given application may not be easy. 

vii. There is a clear need to develop error prediction methodology for the response surface construction 
methods. There is also a need to develop methods to adaptively refine the sampling points based on the 
error estimation. 
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V. Summary 

Response construction methods using Moving Least Squares (MLS), Kriging and Radial Basis Functions (RBF) 
are compared with the Global Least Squares (GLS) method in three numerical examples for derivative generation 
capability. Also, a new Interpolating Moving Least Squares (IMLS) method adopted from the meshless method is 
presented. It is found that the response surface construction methods using the Kriging and RBF interpolation 
yields more accurate results compared with MLS and GLS methods. 
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